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Abstract We study the formation of azimuthons, i.e., rotating spatial solitons, in media 
with nonlocal focusing nonlinearity. We show that whole families of these solutions can be 
found by considering internal modes of classical non-rotating stationary solutions, namely 
vortex solitons. This offers an exhaustive method to identify azimuthons in a given nonlocal 
medium. We demonstrate formation of azimuthons of different vorticities and explain their 
properties by considering the strongly nonlocal limit of accessible solitons. 
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1 Introduction 



There has been growing interest in studies of propagation of optical beams in nonlocal me- 
dia. These are media where the nonlinear response of the material in a specific spatial lo- 
cation is determined not only by the wave intensity in the same location but also in its 
neighborhood. The extent of this neighborhood in comparison to the beam width deter- 
mines the degree of nonlocality. The nonlinear nonlocal response appears to be ubiquitous 
to many physical se ttings. For instance, it is common to media where certain transport pro- 
cesse s such as heat foabbv and Whinnervll 19681 : llitvak etal 19751 : iDavydova and Fishchukl 
Il995h or charge transfer dCalvo et alll2002h. diffusion and/or drift of atomsJSuter and Blas- 
ber dl993nSkupin et al2007l) are responsible for the nonlinearity. It also occurs in systems 
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involvin g long-range interaction of atoms or molecules as it is the case of nematic liquid 
crystals llConti et a 3l2004l|2003h or dipolar Bose-Einstein condensate jGoral et allEoOok 
iNath et all 20071 : Koch et alll2008l) . It has been shown that nonlocal nonlinear response has 
profound consequences on the wave propagation and formation of localized structures (Kro- 
likowski et al 120041) . In particular, nonlocality prevents collapse by providing a stabilization 
mechanis m and enables robust existence of various types of localized s t ructures and spatial 
solitons fcolchugina et ailll980l: iBang et alll2002l ; iBriedis et aj|2005l ; ISkupin et aJ [2006: 
iLashkinl 120071 : lLashkin et alll2007l) . In local nonlinear media the wave perturbation in a 
particular place affects the nonlinearity whic h in turn, influences the wave itself often insti- 
gating its breakup or spatial transformation foesyatnikov et all2005bh . On the other hand, in 
nonlocal media such perturbation is spatially averaged, and hence has a much weaker impact 
on the wave itself, leading to its stabilization. In particular, it has been shown that nonlocal- 
ity support stable propagation of optical vortices, and mul ti-peak solitonic structures wh ich 
are structurally unstable in material with local response dBuccoliero et alll2008L l200~7bl f3). 
A range of particular types of fundamental as well as higher order nonlocal solitons and 
their interactions have bee n even demonstrated e xperimentally in materials with nonlocal 
response of thermal origin dRotschild etaill2006by) . Recently, it has been also shown the- 
oretically that spatial nonlocal response enables realization of the so called azimuthons i.e. 
multiple peak ring-shaped solitons which exhibit angular rotation in propagation (Desyat- 



2005a; lLopez-Aguavo et all l2006bUal: ISkupin et all 12008b . While few specific 
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types of azimuthons have been investigated in various nonlocal model s, using variatio nal 
techniques mentioned above as well as numerical relaxation procedure jLashkinll2008bl l3), 
no general approach to find stable nonlocal azimuthons has been demonstrated so far. In this 
work we study the formation of azimuthons in nonlocal media with gaussian response. We 
show that whole families of these solitons can be tracked down by analyzing bifurcations 
originating from the nonlinear optical potential of vortex solitons. 



2 Azimuthons 

We consider physical systems governed by the two-dimensional nonlocal nonlinear Schro- 
dinger equation 

d 

i—\lf+A ± \if+0Y=O. (1) 
oz 

where d represents the spatially nonlocal nonlinear response of the medium. Its form de- 
pends on the details of a particular physical system. In the following, we will assume that 
the nonlinear response d can be expressed in terms of the nonlocal response function R(r) 

6 = JjR(\r-r'\)\ ¥ (r',z)\ 2 d 2 r' 7 (2) 

where r = xe x +ye y denotes the transverse coordinates. In this work we will use the so-called 
Gaussian model of nonlocality as an illustrative example, 

e = ^//e-^|v/(r',z)|Vr'. (3) 

However, the proposed solutions should exist in many other nonlocal m odels ( iLitvak et all 
ll975l:|Rotschild et alll2006cJ:ISuter and Blasbergll 19931; ISkupin et all|2007 l ; Assanto and Pec 
cianti 



20031: IConti et alll2004|Peccianti et allbOOd : iDenschlag et afeod IFedri and Santosl 
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Azimuthons are a straightforwa rd generalization of the usual ansatz for stationary solu- 
tions (solitons) dPesvatnikov et alll2005ah . They represent spatially rotating structures and 
hence involve an additional parameter, the angular frequency Q (see also ISkrvabin et all 



\ir(r,<t>,z) = U(r,<t>-Qz)e 



(4) 



where U is the complex amplitude function and A the propagation constant. For Q = 0, 
azimuthons become ordinary (nonrotating) solitons. The simplest example of a family of 
azim uthons is the one connecti ng the dipole soliton with the single charged vortex soli- 
ton ( lLopez-Aguavo et~a 3 l2006bh . A single charged vortex consists of two equal-amplitude 
dipole-shaped structures with the relative phase of n/2 representing real and imaginary part 
of U . If these two components differ in amplitudes the resulting structure forms a "rotating 
dipole" azimuthon. If one of the components is zero we deal with the nonrotating dipole 
soliton. In the following we will denote the amplitude ratio of these two vortex components 
by a, which also determines the angular modulation depth of the resulting ring-like structure 
by "1 — a". When higher order (e.g. single charged triple-hump) azimuthons are concerned, 
we can not always identify the angular modulation depth with amplitude ratios of real and 
imaginary part of U. Hence we define the generalized structural parameter a as 



mku \U(r max ,<j>)\ 

CC — 1 

r max i $max ) | 

denotes the coordinates of the maximum value max r ^ \ U 



(5) 



where the tuple (r max ,0 max 

After inserting the ansatz © into Eqs. (QJ and @, multiplying with U* and d^U* resp.. 
and integrating over the transverse coordinates we end up with 



-XM + QL-+I + N = 
-IL 7 +QM' +1' +N' = 0. 



(6a) 
(6b) 



This system relates the propagation constant X and the rotation frequency Q of the az- 
imuthons to integrals over their stationary amplitude profiles, namely 



M 
L z 
I 

M' ■ 
I' = i 



\U(r)\ 2 d 2 r 



I jj [/*(r) — [/(rV 2 r 



U*{r)A ± U(r)d 2 r 

R(\r-r'\)\U{r')\ 2 \U(r)\ 2 d 2 r'd 2 r 
d 2 



d 2 r 



d(j> 



U*(r) 



A ± U(r)d 2 r 



N> = i R(\r-r>\)\U(r')\ 



U(r)d 2 r'd 2 r 



(6c) 
(6d) 
(6e) 
(6f) 

(6g) 
(6h) 
(6i) 



The first two quantities have straightforward physical meanings, namely "mass" (M) and 
"angular momentum" (L z ). We can formally solve for the rotation frequency and obtain (for 
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an alternative derivation see dRozanovll2004h 

_ M(I'+N')-L Z (I + N) 

12 -MM' ■ () 

Note that this expression is undetermined for a vortex beam. For a = 1 [vortex soliton 
V(r) exp(iq(j) + jAqz)], we can assume any value for Q by just shifting the propagation con- 
stant X = ?iq + Q accordingly (Ao accounts for the propagation constant in the non-rotating 
laboratory frame). However, with respect to a particular azimuthon in the limit a — > 1, the 
value of Q is fixed. In what follows, we denote this value by Q \ a =\ ■ 



3 Internal modes and azimuthons 



In this section we will discuss the formation of azimuthons via the process of bifurcation 
from a stationary non-rotating soliton solution, namely a vortex. We assume a certain de- 
formation of the soliton profile while going over from the vortex to azimuthons in the limit 
a — > 1 . Therefore it has to be the shape of vortex deformation which determines Q , since a 
vortex formally allows for all possible rotation frequencies (see the discussion on shifting A 
at the end of Sec.[2]l. 

Let us now look at the azimuthon originating (bifurcating) from a vortex soliton with 
charge q. For this purpose, we recall the eigenvalue problem for internal modes of the non- 
linear potential d which is usually treated in the context of linear stability of nonlinear 
soliton solutions iFirth and Skrvabinll 19971 iDesvatnikov et a 1 2005b). We introduce a small 
perturbation 8V to the vortex soliton V, 



\f/=(V + 8V)e' 



(8) 



plug it into Eqs. (QJ and (|2} and linearize those equations with respect to the perturbation. 
Note that the perturbation 8V(r, (j),z) is complex, whereas the vortex profile V(r) is real 
(w.l.o.g.). The resulting evolution equation for the perturbation 8V is then given by 



d . Id 

«3 Ao + -t— 

dz r dr 



With the ansatz 



dr ) r 2 \d<j) ~^ ^ 



R{\r-r'\)V 2 (r)d 2 r' 8V 
+V II R(\r-r'\)V{r') [8V(r' ,z) + 8V*(r' ,z)]d 2 r' : 



8V = 8Vi { r y ml >>+ iKZ + 8Vi-(r)e- im ^- iK ' z 
we derive the eigenvalue problem for the internal modes 



(9) 



(10) 



1 d 



(m + i 



h+ //fl(|r-r'|)vV)J 2 r' 



8Vi 



+V // R(\r-r'\)V(r') [SVi(r') + SV 2 (r')] cos[m(<j) - <j)')]d 2 r' = kSVi (11a) 



ld_ 

r dr 

-V 



(m-q) 2 



h+ R(\r-r'\)V 2 (r')d 2 r> 



8V 2 



JJ 2?(|r-r'|)V(r') [SV 2 (r') + SVi(r')] cos[m(> - <j>')]d 2 r' = k8V 2 . (lib) 
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Note that since |r — r'| = \Jr 2 + r 12 — 2rr'cos(0 — 0'), all integrals in Jilt are independent 
of <f>. Real-valued eigenvalues of Eq. ( fTTT i (fc = fc*) are termed orbitally stable and the cor- 
responding eigenvector (8V1.8V2) can be chosen as real. If we perturb the vortex V with 
an orbitally stable eigenvector, the resulting wave-function \jf can be written in the form of 
Eq. ((4} with Q = —K/m and A = Ao — qK/m. Thus, it is possible to construct azimuthons in 
the vicinity of the vortex (awl) from 8V: 

U(r,4>)\ z=0 = [V{r) +A r 8V l (r)e im +A r 8V 2 (r)tr im *] e ! '#. (12) 

Used as an initial condition in the propagation equation (0 this object is expected to rotate 
with an angular frequency Q \ a= \ = —K/m. Here, A r > was introduced as the amplitude 
of the perturbation 8V with respect to V. Since we are operating in a linearized system, the 
amplitude of the perturbation as a solution of Eq. flllb is not fixed (just the ratio between 
the components 8V\ and 8V2 is prescribed), but will eventually determine the value of the 
structural parameter a. Generally speaking, the smaller the resulting a the greater the error 
in the constructed initial condition. However, the great robustness of the azimuthons, at least 
in the Gaussian model, allows one to use the initial condition d 12b for quite large perturba- 
tion amplitudes A r . Those strongly perturbed initial conditions result in oscillations of the 
azimuthon upon propagation. However, the azimuthon is structurally stable and does not 
decay into other soliton solutions like the single-hump ground state. Moreover, such initial 
conditions play a role of excellent "initial guesses" for solver routines to find numerically 
exact azimuthons. 

4 Higher order azimuthons 

In a recent publication we have used the approach presented above to characterize the ro- 
tating dipole a zimuthon, which c onnects the single charged vortex (q = 1) to the stationary 
dipole soliton dSkupin et alll2008l) . As already mentioned there, solving the eigenvalue prob- 
lem dl 11 1 can be used as an exhaustive method for finding families of azimuthons which 
originate from a vortex soliton. However, it should be stressed that not all orbitally sta- 
ble eigenvalues can be linked to a family of azimuthons. This is obvious for eigenval- 
ues with |fc| > Ao in the continuous part of the spectrum. Hence, we can conclude that 
|X2 1 cr=i I < \Xo\/m for azimuthons in the vicinity of the vortex (a ~ 1). Note that the param- 
eter m determines the number of humps of the rotating structure. 

Looking at a higher order azimuthon, e.g., a single charged rotating triple hump (q = 
1, m = 3), the natural question one may pose is whether this particular family of azimuthons 
is connect ed to a second bifurcation from a vortex soliton, as pr edicted by variational cal- 
culations dDesvatnikov et "a3 l2005al : II opez-Aguavo et "all l2006ah . In the case of the rotat- 
ing double hump we had a natural candidate, namely the stationary dipole; the existence 
of an analogous solution like a stationary tripole is not evident. It turns out that the rotat- 
ing triple hump azimuthon with lowest absolute rotation frequency Q connects single and 
double charged vortex with opposite sign of charge, in contrast to variational predictions 
dLopez-Aguavo et "all l2006ah . At least in the highly nonlocal regime the rotating frequency 
does not change much when we follow the family with constant mass (here M = 630 and 
— 3 > Q > —3.5), we do not find a solution with Q = 0. What we do find is a solution with 
vanishing angular momentum L z , because the two limiting vortices have opposite sign of 
charge. 

Analysis of internal modes of both vortices (V) (charge q = 1 and q = —2, charge of 
perturbation m = 3) reveals eigenvalues and eigenvectors where the family of azimuthons 
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Fig. 1 Internal modes (m = 3) of the single charges vortex (q = 1) with eigenvalue K = 10.3 (left) and the 
double charged vortex (q = —2) with eigenvalue K = 9.1 (right). SVo clearly dominates the nearly invisible 
8V{ component in the left picture, whereas the opposite is the case in the right one. 
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Fig. 2 (First row) Tracking the family starting from the single charged vortex (q = 1), the respective angular 
momenta and structural parameter are (from left to right) L- = 450, a = 0.65; L z = 270, a = 0.53; and L z = 0, 
a = 0.4. (Second row) Tracking the family starting from the double charged vortex (q = —2); L z = — 1 100, 
a = 0.55; L z = -770, a = 0.3; and L z = 0, a = 0.4 (from left to right). 



emerges. Figure Q] shows the two eigenvectors for mass M = 630, the resulting rotation fre- 
quency is Q = — k/3. Now we can track the family starting from both vortices (q = 1 and 
q = —2) perturbed with the appropriate eigenvectors and constant mass M = 630. Equa- 
tion J 12b serves as initial condition (we increase the perturbation amplitude) to a Newton 
solver. We follow both branches till we reach L z = 0, where the two solutions coincide (see 
last plot in each row of Fig. [2j». Interestingly, we observe that the triple hump azimuthon 
with L z = is not the one with maximum modulation depth I — a. 

From a topological point of view, the above findings are somehow surprising because 
the two limiting vortices have different charge (q = 1 and q = —2). It is possible to under- 
stand this interesting feature when looking at the azimuthon close to the respective vortices. 
Starting from the double charged vortex (q = —2), the azimuthon is created by adding a 
counter-rotating single charged vortex (m + q = 1, see right panel in Fig.|TJ. However small 
the amplitude of this single charged vortex might be, in the vicinity of the origin it will al- 
ways be dominant as it grows as ~ r, whereas the double charged vortex as ~ r 2 . Thus, the 
azimuthon has aq = I vorticity at the origin, and on a ring where the amplitudes of the two 
vortices are equal lie three phase singularities with charge — 1 (see Fig. [3] for a schematic 
sketch). As we can see in Fig. [2] the radius of this ring grows when we follow the family of 
azimuthons towards the single charged vortex, and the three singularities with total charge 
—3 move far away from the origin and finally disappear when we approach the vortex with 
total charge q = 1 (see discussion below). It is important to note that these three phase sin- 
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Fig. 3 Sketch of the co-rotating topological charges observed in the azimuthon family of Fig. [2] 



gularities have fixed positions with respect to the position of the three humps and follow the 
amplitude rotation of the azimuthon (co-rotating). 

We will now discuss in a greater detail the three co-rotating phase singularities, in partic- 
ular, how they disappear when we approach the vortex with total charge q = 1. To this end, 
we analyze the asymptotic behavior for large r of the three components of the azimuthon V, 
8V\ and 8V2 [see Eq. l[T2"lll. For sufficiently large r, the convolution term 8 in Eqs. Q and 
dTTT) can be neglected when compared to the term ~ 1 /r 2 in the transverse Laplacian. Then, 
using the modified Bessel functions, one can find the asymptotic behavior of the involved 
functions easily: 




l 4(m + qf-\ 
4(m-q) 2 - 1 



1 + „ == + e\ 



(13a) 
(13b) 
(13c) 



To find the radius r s where the phase singularities appear, one has to equal the amplitudes in 
the following manner: 

\V(r)\=A r \SVi(r) + SV 2 (r)\. (14) 

It is obvious from Eqs. dl3ab that such a radius exists for arbitrary small A r , because one of 
the 5V, decays always slower than V for r — > 00. In our example we have K > 0, and therefore 
8V2 is responsible for creating our three co-rotating phase singularities. For A r — > we 
find that r s — » 00, the singularities move to infinite distances form the origin and (formally) 
vanish for A r = 0. However, for practical observations in, e.g., numerical simulations those 
co-rotating phase singularities become irrelevant when the surrounding amplitude becomes 
small. 

The observation that the above triple hump azimuthon has almost constant angular fre- 
quency when we follow the family for constant mass M regime can be explained going 
over to the highly nonlocal limit. This rotation is not a purely nonlinear phenomenon, but 
is mainly a consequence of mode beating. Let us have a look at the related linear limit 
where we replace the non local response © by the Ga ussian kernel times mass M (similar 
to Snyder-Mitchel model dSnvd er and Mitc hell 1997)). In this linear problem we can find 
several eigenmodes (see Fig.|4]l. Mode beating between single {q = 1) and double {q = —2) 
charged vortices predicts Q ~ (A C4 — Xq 2 )/3 = —4.2, which is not too far from the rotation 
frequency observed in the nonlocal nonlinear problem. 

Another evidence that the linear contribution to the rotation of the triple hump dominates 
is that Q increases strongly with mass M (and A), as expected from Q ~ (Ac 4 — Ac 2 ) /3. 
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Fig. 4 First six linear modes of a Gaussian potential with mass 630: Ground state G\ (r), Xg t = 86.6; sin- 
gle charged vortex Gi{r) exp(i</>), Xq = 73.5; humped ring state Gi(r), Xq^ = 61.5; double charged vor- 
tex G4(r)exp(i20), Aq 4 = 60.9; single charged double ring G5(r)exp(i0), Xq 5 = 50; triple charged vortex 
G(,(r) exp(f'30), Xq 6 = 48.9 (from left to right, top to down). 



E.g., for M = 200 we find Q. ~ —1.2. In contrast to that, the double hump azimuthon con- 
necting singl e charged vortex a nd stationary dipole shows almost no dependency of Q 
on the mass dSkupin et al l2008). This resembles the fact that in the linear problem men- 
tioned above, mode beating pr edicts Q ~ (Xq-, — Xc n )l 2 = for this structure. In fact, fol- 
lowing a reasoning similar to iBuccoliero et all <J2008), we can identify in the expression 
I0 for the rotation frequency a linear and nonlinear contribution, Q = £2i„ +Q n i n . In the 
limit © = Mexp(— r 2 /2)/2n and U a superposition of "linear" modes we readily see that 
MN' — L.N = 0, and any rotation is due to 

MI' - LJ 

n >- = ij^mm- ^ 

In the special case that U is a superposition of degenerated linear modes, we find MI' — LJ = 
and thus Q/„ = 0. However, if we consider the original nonlinear system where © is given 
by Eq. ©, an additional (nonzero) nonlinear contribution 

n nln = MN '- L ^ N (16) 

L\ — MM 

to the rotation frequency occurs. 

Once we have computed the internal mode of a vortex, we can construct all azimuthons 
branching from it. For example, Fig. [5] shows a rotating five-hump azimuthon emanating 
from our double charged vortex (q = — 2). The corresponding internal mode shows a typical 
r ( 5 - 2 ) dependence near the origin in 8V\ , the amplitude of 8V2 is very small. Hence, for the 
azimuthon, we see the double charged phase singularity {q = —2) of the vortex in the origin. 
As observed for the rotatin triple-hump above, five singularities with q = 1 appear on a ring, 
and they are expected to move inwards when we follow the azimuthon family towards the 
triple charged vortex with q = 3. 

We can also easily identify families of azimuthons previously found using a special 
ansatz. For instance, the double charged vortex (q = —2) shows another internal mode for 
m = 2 with eigenvalue K = —3.1. As can be seen in Fig.[6l the resulting azimuthon belongs to 
the family connecting Hermite-Gaussian and Laguerre-Gaussian self-trapped modes HN20 
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Amplitude 



Fig. 5 Internal mode (in = 5) of the double charged vortex (q = —2) with eigenvalue K = —8.6, and corre- 
sponding emanating azimuthon (M = 630, a = 0.56). SVi clearly dominates the small 8V2 component in the 
left picture. 




Fig. 7 Internal mode (m = 2) of the groundstate (q = 0) with eigenvalue K = —23.3, and corresponding 
emanating azimuthon (M = 630, £2 = 14.3). 



and LN20 teuccoliero et all [2008b . Note that for this solution our definition of the struc- 
tural parameter a does not make sense, hence in the caption of Fig. [6] we give the rotation 
frequency Q instead, to characterize the azimuthon. 

Last but not least, we want to note here that the concept of azimuthons branching from 
solitons is not limited to vortices. E.g., the single-hump ground state (M = 630, q = 0) fea- 
tures am = 2 internal mode with K = — 23.3. The emanating azimuthon looks like a rotating 
bone, and possesses two phase singularities with opposite charges on an axis perpendicular 
to the "bone" axis (see Fig. [7}. The very high rotation frequency is again a manifestation of 
linear mode beating, we find Q ~ (Ag, — Ag 4 ) /2 = 12.9. 

If we reduce mass and therefore, in our scaling, reduce nonlocality these results may 
change. First of all, solitons and azimuthons are expected to become unstable. Moreover, 
we observe that the second component of the perturbation becomes larger in amplitude 
when we leave the highly nonlocal regime (note that it is almost invisible in Figs.[T]|5]|5] 
and |7]i. Also, certain types of internal modes may vanish or new ones may appear. Hence, 
careful analysis of the internal spectra of soliton solutions is necessary to predict structure 
of azimuthons in a given regime or model, e.g., the local nonlinear Schrodinger equation. 
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5 Conclusion 

In conclusion, we have demonstrated a simple method for identifying rotating solutions in 
nonlocal nonlinear media. We computed azimuthon solutions and their rotation frequen- 
cies numerically and showed that in the limit of minimal azimuthal amplitude modulation, 
i.e., close to a vortex soliton, the rotation frequency is determined uniquely by eigenvalues 
of the bound modes of the linearized version of the respective stationary nonlocal solu- 
tion. Moreover, the intensity profile of the resulting azimuthons can be constructed from the 
corresponding linear eigensolution. This offers a straightforward and exhaustive method to 
identify rotating soliton solutions in a given nonlinear medium. At least, in the highly non- 
local regime, we find families of azimuthons which connect vortex solitons with different 
topological charge. 
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